function indicator = featureSelection(featureCalculator, dim)
% Visual feature selection (d to d', d' < d) via trace ratio criterion.
% Input:
%   featureCalculator: a FeatureCalculator instance, the features computed
%     by which is a n-by-d matrix;
%   (coordinates: n-by-p matrix of ground-truth coordinates)
%   dim: number of selected dimensions d'
%
% Output:
%   indicator: 1-by-d vector, where 1 indicates a selected feature.

	tau = 80;  % the threshold of pose distance, in millimetre(mm)
	sigma2 = 5 * 5;  % parameter sigma = 5
	nPairs = 10000;
	nData = featureCalculator.imageCollection.nImages;
	% The number of columns is nData - 1, in order to make sure the two
	% fields of any pair wouldn't be the same.
	% (1) [4] [7]    [4] [7]
	% [2] (5) [8] => [2] [8]
	% [3] [6] (9)    [3] [6]
	nElementsOfPairs = nData * (nData - 1);
	idx = sort(randsample(nElementsOfPairs, nPairs));
	[p1 p2] = ind2sub([nData, nData - 1], idx);  % pairs
	p2 = p2 + (p1 <= p2);  % moves to right if above the diagonal
	
	% The following line is equivalent to:
	%   f1 = features(p1, :);
	%   f2 = features(p2, :);
	%   ps1 = coordinates(p1, :);
	%   ps2 = coordinates(p2, :);
	% where `features` is the whole features that featureCalculator computes, 
	% and `coordinates` is the corresponding coordinates.
	[f1, f2, ps1, ps2] = featureCalculator.calculatePair(p1, p2, true);

	df = f1 - f2;
	dp = poseDistance(ps1, ps2);  % nPairs-by-1


	% ----8<--------------------
	dpMax = max(dp);

	pairwiseSimilar = exp(-dp / sigma2);

	pairwiseDissimilar = dp;
	ind0 = pairwiseDissimilar < tau;
	pairwiseDissimilar(ind0) = 0;
	pairwiseDissimilar(~ind0) = (pairwiseDissimilar(~ind0) / dpMax) .^ 2;
	
	% -------------------->8----


	% The following code block vectorizes the commented one for speed:
	% Uw = 0; Ub = 0;
	% for i = 1:nPairs
	% 	d_ = df(i, :);
	% 	D = d_' * d_;
	% 	Uw = Uw + D * pairwiseSimilar(i);  % within
	% 	Ub = Ub + D * pairwiseDissimilar(i);  % between
	% end
	
	% (av)' * (av) = (a^2) (v' * v)
	
	pairwiseSimilar = pairwiseSimilar .^ 0.5;  % squared in multiplication
    pairwiseDissimilar = pairwiseDissimilar .^ 0.5;
	dw = bsxfun(@times, df, pairwiseSimilar);
	db = bsxfun(@times, df, pairwiseDissimilar);
	Uw = dw' * dw;
	Ub = db' * db;
	
	indicator(1, size(df, 2)) = false;  % init
	indices = graphScore(Ub, Uw, dim);
	indices = indices(1:dim);  % A bug in `graphScore`?
	indicator(indices) = true;
end


